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Cosmic rays are atomic nuclei with velocities close to the speed 
of light. Ultra-high energy cosmic rays (UHECRs) are approximately 
ten million times more energetic than the most extreme particles 
produced at the Large Hadron Collider and likely originate from 
nearby extragalactic sources. Important astrophysical questions in- 
clude: what phenomenon accelerates particles to such large energies, 
which astronomical objects host the accelerators, and what sorts of 
nuclei end up being energized? We develop a multilevel Bayesian 
framework for assessing evidence for association of UHECRs and can- 
didate source populations that has four levels: (1) a source property 
description level, including specification of a candidate source popu- 
lation and the distribution of cosmic ray intensities among sources; 
(2) a cosmic ray production level, using a marked Poisson point pro- 
cess to model latent properties (arrival times and energies) of cosmic 
rays; (3) a propagation level, modeling deflection of cosmic ray tra- 
jectories by cosmic magnetic fields; (4) a detection and measurement 
level, accounting for detector efficiency and exposure, and measure- 
ment errors in arrival directions and energies. Our framework can 
flexibly accommodate astrophysical components anywhere on a spec- 
trum from simplicity to great complexity and realism. The number of 
possible associations is huge even for modest-sized datasets. We de- 
scribe a Markov chain Monte Carlo algorithm for estimating model 
parameters and comparing models by computing, via Chib's method, 
marginal likelihoods and Bayes factors. 

We demonstrate the framework by using simple model components 
to analyze observations of 69 UHECRs by the Pierre Auger Obser- 
vatory (PAO) during the years 2004-2009, using a volume-complete 
catalog of 17 nearby active galactic nuclei as a candidate host popu- 
lation. The reported data are incomplete; an early portion of the data 
("Period 1," with 14 events) was used to set an energy cut maximizing 
a measure of anisotropy in that data; only data for cosmic rays with 
energies above that cut are reported. Also, measurement errors are 
approximately summarized. These factors are problematic for inde- 
pendent analyses of PAO data. Within the context of "standard can- 
dle" source models (all sources with the same isotropic cosmic ray 
emission rate), and considering only the untuned data subsequent 
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to Period 1, there is no significant evidence favoring association of 
UHECRs witli nearby active galactic nuclei (AGN) vs. attributing 
them to sources drawn from an isotropic background distribution. 
The highest-probability associations are with the two nearest AGN, 
Centaurus A and its neighbor, NGC 4945. If the association model 
is adopted, the fraction of UHECRs that may be associated is likely 
nonzero but is constrained to be well below 50%. Relatively modest 
magnetic deflection angular scales of « 3° to 30° are favored. Mod- 
els that assign a large fraction of UHECRs to a single nearby source 
(e.g., Centaurus A) are ruled out unless very large deflection scales 
are specified a priori, and even then they are disfavored. However, 
including the Period 1 data alters the conclusions significantly, and a 
simulation study supports the idea that the Period 1 data are anoma- 
lous, presumably due to the tuning. Accurate and optimal analysis of 
future data will likely require more complete disclosure of the data. 

1. Introduction. Cosmic ray particles are naturally produced, posi- 
tively charged atomic nuclei arriving from outer space with velocities close 
to the speed of light. Their celestial origin was discovered at the beginning of 
the 20th century using electrometers hoisted by balloons into the upper at- 
mosphere [22]. Nearly a century's subsequent observations have determined 
that the impinging flux includes particles with energies spanning a dozen or 
more orders of magnitude [15]. 

The origin of cosmic rays is not well-understood. The Lorentz force expe- 
rienced by a charged particle in a magnetic field alters its trajectory. Simple 
estimates imply that cosmic rays with energy E < 10^^ eV have trajectories 
so strongly bent by the Galactic magnetic field that they are largely trapped 
within the Galaxy.^ The acceleration sites and the source populations are 
not definitively known but probably include supernovae, pulsars, stars with 
strong winds, and stellar-mass black holes. For recent reviews, see [15, 23]. 
More mysterious, however, are the highest energy cosmic rays. 

In the 1960s, large arrays of cosmic ray detectors began detecting cosmic 
rays with enormous energies; by 1991 a few events were seen with energies 
~ 100 EeV (where 1 EeV = 10^® eV). In the 1990's cosmic ray astronomers 
built huge detector arrays specifically targeting these energetic events. The 
largest and most recent of these were the Akeno Giant Air Shower Ar- 
ray (AGASA; [13]) and the High Resolution Fly's Eye (HiRes; [10]); each 
detected a few dozen cosmic rays with E > 10 EeV, so-called ultra-high 
energy cosmic rays (UHECRs). For recent reviews, see [26, 27]. The sources 
of these cosmic rays likely form a relatively nearby extragalactic population. 

^We follow the standard astronomical convention of using "Galaxy" and "Galactic" 
(capitalized) to refer to the Milky Way galaxy. 
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On the one hand their trajectories are only weakly deflected by galactic 
magnetic fields so they are unconfined to the galaxy from which they orig- 
inate. On the other hand, they are unlikely to reach us from distant (and 
thus isotropically distributed) cosmological sources. Cosmic rays with ener- 
gies above the Greisen-Zatsepin-Kuzmin (GZK) limit of ~ 50 EeV should 
scatter off of cosmic microwave background photons, losing some of their 
energy to pion production with each interaction [19, 40]. Thus the universe 
is not transparent to UHECRs; they are not expected to travel more than 
about 50 to 100 megaparsecs (Mpc) before their energies fall below the GZK 
limit. Detectable UHECRs likely emanate from relatively nearby extragalac- 
tic sources. Notably, over this cosmologically modest distance scale there is 
significant anisotropy in the distribution of matter that should be reflected 
in the arrival directions of UHECRs (albeit distorted by magnetic deflec- 
tion) . Astronomers hope that continued study of the directions and energies 
of UHECRs will address the fundamental questions of the field: What phe- 
nomenon accelerates particles to such large energies? Which astronomical 
objects host the accelerators? What sorts of nuclei end up being energized? 
In addition, UHECRs probe galactic and intergalactic magnetic fields. 

The flux of UHECRs is very small, approximately 1 per square kilometer 
per century for energies E>50 EeV. Large detectors are needed to find these 
elusive objects. The construction of the largest and most sensitive detector 
to date, the Pierre Auger Observatory (PAO) [4] in Argentina, has marked a 
significant step towards that goal. The observatory uses a combination of air 
fluorescence telescopes and water Cerenkov surface detectors to observe the 
air shower generated when a cosmic ray interacts with nuclei in the upper 
atmosphere over the observatory. The surface detectors (SDs) operate con- 
tinuously, detecting energetic subatomic particles produced in the air shower 
and reaching the ground. The fluorescence detectors (FDs) image light from 
the air shower and supplement the surface detector data for events detected 
on clear, dark nights.^ PAO began taking data in 2004 even as it was under 
construction; it reached its baseline design in June 2008 with an array of 
~ 1600 SDs covering ~ 3000 km^, surrounded by four fluorescence telescope 
stations (with six telescopes in each station) observing the atmosphere over 
the array. 

By 31 August 2007, PAO had detected 81 UHECRs with E > AO EeV (see 

^The FD on-time is about 13% [3], but analysis can reveal complications preventing 
use of the data — e.g., obscuration due to light cloud cover, or showers with significant 
development underground — so fewer than 13% of events have usable FD data. These few 
so-called hybrid events are important for calibrating energy measurements and provide 
information about cosmic ray composition vs. energy. 
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[31], hereafter PAO-07), finding clear evidence of an energy cutoff resembling 
the predicted GZK cutoff, i.e., a sharp drop in the energy spectrum above 
~ 100 EeV and a discernable pile-up of events at energies below that [3]. This 
supports the idea that the UHECRs must originate in the nearby universe. 

The PAO team searched for correlations between the cosmic ray arrival 
directions and the directions to nearby active galactic nuclei (AGN) (initial 
results were reported in PAO-07; further details and a catalog of the events 
are in [32], hereafter PAO-08). AGN are unusually bright cores of galaxies; 
there is strong (but indirect) evidence that they contain supermassive black 
holes rapidly accreting nearby gas and stars and ejecting some of the accret- 
ing material in energetic, jet-like outflows. AGN are theoretically favored 
sites for producing UHECRs; electromagnetic observations indicate particles 
are accelerated to high energies near AGN. The PAO team's analysis was 
based on a significance test that counted the number of UHECRs with best- 
fit directions within a critical angle, tp, of an AGN in a catalog of local AGN 
candidate sources (more details about the catalog appear below); the num- 
ber was compared with what would be expected from an isotropic UHECR 
directional distribution using a p-value. A simple sequential approach was 
adopted. The earliest half of the data was used to tune three parameters 
defining the test statistic by minimizing the p-value. The parameters were: 
the critical angle ip; a maximum distance cutoff, -Dmax, for possible hosts; 
and an energy threshold, E'thi with only UHECRs with E > E'th considered 
to be associated with AGN. With these parameters tuned (£'th = 56 EeV, 
i; = 3.1°, D 

max — 75 Mpc), the test was applied to the later half of the 
data; 13 UHECRs in that period had E > £^th- The resulting p- value of 
1.7 X 10"'^ was taken as indicating the data justify rejecting the hypothesis 
of isotropic arrival directions "with at least a 99% confidence level." The 
PAO team was careful to note that this result did not necessarily imply that 
UHECRs were associated with the cataloged AGN, but rather that they 
were likely to be associated with some nearby extragalactic population with 
similar anisotropy. 

Along with these results, the PAO team published a catalog of energy and 
direction estimates for the 27 UHECRs satisfying the E > ii^th criterion, in- 
cluding both the earliest 14 events used to define E'th, and the 13 subsequent 
events used to obtain the reported p-value (the PAO data are proprietary; 
measurements of the other 54 events used in the analysis were not pub- 
lished). Their exciting but admittedly suggestive statistical result spurred 
subsequent analyses of these early published PAO UHECR arrival direc- 
tions, adopting different methods and aiming to make more specific claims 
about the hosts of the UHECRs. Roughly speaking, these analyses found 
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similarly suggestive evidence for anisotropy, but no conclusive evidence for 
any specific association hypothesis. 

In late 2010, the PAO team published a revised catalog, including new 
data collected through 2009 ([33]; hereafter PAO-10). An improved analysis 
pipeline revised the energies of earlier events downward by 1 EeV; accord- 
ingly, the team adopted £'th = 55 EeV on the new energy scale. The new 
catalog includes measurements of 42 additional UHECRs (with E > -Eth) 
detected from 1 September 2009 through 31 December 2010, for a total of 69 
events. A repeat of the previous analysis (adding the new events but again 
excluding the early tuning events) produced a larger p- value of 3 x 10"^, 
i.e., weaker evidence against the isotropic hypothesis. The team performed 
a number of other analyses (including considering new candidate host pop- 
ulations). Despite the growth of the post-tuning sample size from 14 to 55, 
they found evidence for anisotropy weakened, across diverse analyses. Time- 
resolved measures of anisotropy provided puzzling indications that later data 
might have different directional properties than early data, although the 
sample size is too small to demonstrate this conclusively. 

Here we describe a new framework for modeling UHECR data based on 
Bayesian multilevel modeling of cosmic ray emission, propagation, and de- 
tection. Astrophysically, a virtue of this approach is that physical and exper- 
imental processes have explicit representations in the framework, facilitat- 
ing exploration of various scientific hypotheses, and physical interpretation 
of the results. This is in contrast to hypothesis testing approaches where 
elements such as angular and energy thresholds only implicitly represent 
underlying physics, and potentially conflate astrophysical and experimental 
effects (e.g., magnetic scattering of trajectories, and measurement errors in 
direction). Statistically, an important virtue of our framework is the abil- 
ity to handle a priori uncertainty in model parameters via marginalization. 
Marginalization also accounts for the uncertainty in such parameters via 
weighted averaging, rather than fixing them at precise, tuned values. This 
eliminates the need to tune energy, angle, and distance scales with a subset 
of the data that must then be excluded from a final analysis. Such param- 
eters are allowed to adapt to the data, but the "Ockham's razor" effect 
associated with marginalization penalizes models for fine-tuned degrees of 
freedom, thereby accounting for the adaptation. 

In this paper we describe our general framework, computational algo- 
rithms for its implementation, and results from analyses based on a few 
representative models. Our models are somewhat simplistic astrophysically, 
although similar to models adopted in previous studies. We do not aim to 
reach final conclusions about the sources of UHECRs; the focus here is on 
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developing new methodology and demonstrating the capabilities of the ap- 
proach in the context of simple models. 

An important finding is that thorough and accurate independent analy- 
sis of the PAO data likely requires more data than has so far been publicly 
released by the PAO collaboration. In particular, although our Bayesian 
approach eliminates the need for tuning, in the absence of publicly avail- 
able "untuned" data (i.e., measurements of lower-energy cosmic rays), we 
cannot completely eliminate the effects of tuning from analyses of the pub- 
lished data (Bayesian or otherwise). Additionally, a Bayesian analysis can 
(and should) use event-by-event (i.e., heteroskedastic) measurement uncer- 
tainties, but these are not publicly available. Finally, astrophysically plau- 
sible conclusions about the sources of UHECRs will require models more 
sophisticated than those we explore here (and those explored in other re- 
cent studies). The PAO observatory continues to operate; we hope future 
PAO publications will include more complete event catalogs, enabling more 
thoroughgoing analyses than are presently possible. 

2. Description of cosmic ray and candidate host data. 

2.1. Cosmic ray data. The reported PAO measurements depend not only 
on the intrinsic particle population but also on many experimental and al- 
gorithmic choices in the detection and analysis chain, many of them associ- 
ated with the need to distinguish between events of interest and background 
events from uninteresting but uncontrollable sources (e.g., natural radioac- 
tivity). UHECRs can impinge on the observatory at any time, from any 
direction and with any energy. However, virtually no background sources 
produce events with properties mimicking those of very high energy cos- 
mic rays arriving from directions well above the horizon. Cosmic rays with 

> 3 EeV arriving from any direction lying within a large window on the 
sky create air showers detected with nearly 100% efficiency (no false posi- 
tives, no false dismissals). The SDs and FDs measure the spatio-temporal 
development of the air shower which allows the energy and arrival direction 
to be measured. The uncertainties depend upon how many counters of each 
type are triggered plus the systematic and statistical uncertainties implicit 
in modeling the development of the air shower. The PAO team reports en- 
ergy and arrival direction estimates for each cosmic ray falling within the 
geometric bounds of its zone of secure detection.^ 

^The directional criterion adopted for the PAO catalogs is that an event is reported if 
it best-fit arrival direction is within 60° of the observatory's zenith, the local normal to 
Earth's surface at the time of the event. 
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We consider the Nq = 69 UHECRs with energies E > Eth = 55 EeV 
cataloged in PAO-10, which reports measurements of ah UHECRs seen by 
PAO through 31 December 2009 with E > £^th) based on analysis of the 
surface detector data only. Although our framework does not tune an event 
selection criterion, for interpreting the results it is important to remember 
that the ii^th = 55 EeV threshold value was set to maximize a signature 
of anisotropy in an early subset of the data. The tuning data included the 
14 earliest reported events, detected from 1 January 2004 to 26 May 2006 
(inclusive; Period 1), as well as numerous unreported events with E < E^h- 
The first published catalog in PAO-08 included 13 subsequent UHECRs ob- 
served through 31 August 2007 (Period 2). The PAO-10 catalog includes 42 
additional UHECRs observed through 31 December 2009 (Period 3). Table 1 
in PAO-10 provides information about the three periods, including the sky 
exposure for each period, which is not simply proportional to duration (the 
observatory grew in size considerably through 2008). Data for cosmic rays 
with E < Eth are not publicly available.^ 

The direction estimate for a particular cosmic ray is the result of a compli- 
cated analysis of time series data from the array of PAO surface detectors.^ 
Roughly speaking, the direction is inferred by triangulation. The analysis 
produces a likelihood function for the cosmic ray arrival direction, lo (a unit 
vector on the celestial sphere) . The shapes of the likelihood contours are not 
simple, but they are roughly azimuthally symmetric about the best-fit direc- 
tion. The PAO-10 catalog summarizes the likelihood function with a best-fit 
direction, and a typical directional uncertainty of ~ 0.9° corresponding to 
the angular radius of an azimuthally symmetric 68.3% confidence region. We 
use these summaries to approximate the likelihood functions with a Fisher 
distribution with mode at the best-fit direction for each cosmic ray, and 
with concentration parameter Kc = 9323, corresponding to a 68.3% confi- 
dence region with an angular radius of 0.9°. Let dj denote the data associated 
with cosmic ray i, and Wj denote its actual arrival direction. The likelihood 
function for the direction is 

(1) i{0Ji) ■= P{di\uii) W r exp^KcUi ■ Ui), 

Att smh(Kcj 

where rij denotes the best-fit direction for cosmic ray i, and we have scaled 
the likelihood function so its integral over ojj is unity, merely as a convenient 

*The PAO web site hosts public data for 1% of lower-energy cosmic rays, but the sample 
is not statistically characterized and UHECRs are not included. 

^The SD data may be supplemented by data from the fluorescence detectors for hybrid 
events observed under favorable conditions, but there are very few such events at ultra-high 
energies, and PAO-10 reports analysis of SD data only. 
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convention. Bonifazi et al. [9] provide more information about the PAO di- 
rection measurement capability. Note that the expected angular scale of 
magnetic deflection is larger than the PAO directional uncertainties, signif- 
icantly so if UHECRs are heavy nuclei (see § 3.4). 

Similarly, the analysis pipeline produces energy estimates for each event. 
These estimates have significant random and systematic uncertainties. The 
PAO-10 catalog reports a best-fit energy for each reported event, with a 
typical energy resolution of ~ 15% [3]. This uncertainty is the random un- 
certainty in estimates based on surface detector data.^ The energy scale of 
UHECRs is so far beyond the scales of easily measured natural phenomena 
that an accurate absolute energy calibration for UHECRs requires signifi- 
cant extrapolation of particle physics theory to regimes not explored directly 
by experiment, and this is a significant source of systematic (correlated) un- 
certainty in the measured energies. A key feature of PAO is the presence of 
both surface and fluorescence detectors at the same site; this enables a more 
secure calibration than was possible for earlier experiments that relied on 
a single type of detector. Analysis and modeling of FD data produces en- 
ergy estimates with an absolute energy scale systematic uncertainty of 22%. 
The absolute energy scale for the SD analysis is calibrated by comparison 
of the SD and FD estimates for hybrid events, and the resulting systematic 
uncertainty in the SD energy scale with respect to the FD scale is ~ 10%. 

The models we study here do not make use of the reported energies and 
are unaffected by these uncertainties. But our framework readily general- 
izes to account for energy dependence. In principle it is straightforward to 
account for the random uncertainties, but a consistent treatment requires 
data for events below any imposed threshold: the true energies of events 
with best-fit energies below threshold could be above threshold (and vice 
versa for those with best-fit energies above threshold); accounting for this 
requires data to energies below astrophysically important thresholds. The 
systematic uncertainties become important for joint analyses of PAO data 
with data from other experiments, and for linking results of spectral analyses 
to particle physics theory. 

2.2. Candidate source catalog. As candidate sources for the PAO UHE- 
CRs, the analysis reported in PAO-07 and PAO-08, and several subsequent 
analyses, considered 694 AGN within ~ 75 Mpc from the 12th catalog assem- 
bled by Veron-Cetty and Veron [37] (VCV). This catalog includes data on all 

®For hybrid events, the uncertainty can be better than 6%. The percentages indicate 
the sizes of the root-mean-square errors of energy estimates in caUbration experiments; the 
errors appear approximately normally distributed and have some dependence on cosmic 
ray composition [1]. 
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AGN and quasars (AGN with star-like images) with pubhshed spectroscopic 
redshifts; it includes observations from numerous investigators using diverse 
equipment and AGN selection methods, and does not represent a statisti- 
cally well-characterized sample of AGN7 Subsequent analyses in PAO-10, 
and a few other analyses, used more recent catalogs of active galaxies or 
normal galaxies, including flux-limited catalogs (i.e., well-characterized cat- 
alogs that contain all bright sources within a specified volume, but dimmer 
sources only in progressively smaller volumes). 

For the representative analyses reported here, we consider the 17 AGN 
cataloged by Goulding et al.[18] (2010; hereafter GIO) as candidate sources. 
This is a well-characterized ?;o/Mme-limited sample; it includes all infrared- 
bright AGN within 15 Mpc. For each AGN in the calalog, we take its po- 
sition on the sky, Wk (A; = 1 to Ns), and its distance, D^, to be known 
precisely (galaxy directions have negligible uncertainties compared to cos- 
mic ray directions). Notably, this catalog includes Centaurus A (Cen A), 
the nearest AGN [D ~ 4.0 Mpc), an unusually active and morphologically 
peculiar AGN. Theorists have hypothesized Cen A to be a source of many 
or even most UHECRs if UHECRs are heavy nuclei, which would be de- 
flected through large angles; see [7, 17, 8]. The small size of this catalog 
facilitates thorough exploration of our methodology: Markov chain Monte 
Carlo algorithms can be validated against more straightforward algorithms 
that could not be deployed on large catalogs, and simulation studies are fea- 
sible that would be too computationally expensive with large catalogs. Also, 
for simple "standard candle" models (adopted here and in other studies), 
that assign all sources the same cosmic ray intensity, little is gained by con- 
sidering large catalogs, because assigning detectable cosmic ray intensities 
to distant sources would imply cosmic ray fluxes from nearby sources too 
large to be compatible with the data. 

We also include an isotropic background component as a "zeroth" source. 
This allows a model to assign some UHECRs to sources not included in the 
AGN catalog (either galaxies not cataloged, or other, unobserved sources). 
In addition, we consider an isotropic source distribution for all cosmic rays 
(i.e., a model with only the zeroth source) as a "null" model for compari- 
son with models that associate some cosmic rays with AGN or other discrete 
sources. An isotropic distribution is convenient for calculations and has been 
adopted as a null hypothesis in several previous studies. Historically, before 
PAO's convincing observation of a GZK-like cutoff in the UHECR energy 

^VCV say of the catalog, "This catalogue should not be used for any statistical analysis 
as it is not complete in any sense, except that it is, we hope, a complete survey of the 
literature." 
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spectrum, the isotropic distribution was meant to represent a distant cosmo- 
logical origin for UHECRs. Accepting the null would indicate that the GZK 
prediction was incorrect, and that changes in fundamental physics would be 
required to explain UHECRs. In light of PAO's compelling observation of a 
GZK-like cutoff (with its implied ~ 100 Mpc distance scale), interpreting an 
isotropic null or background component is problematical if there are many 
light nuclei among the UHECRs. We adopt it here both for convenience and 
due to precedent. We discuss this further below. 

2.3. Sky maps. Figures 1 and 2 show sky maps displaying the directions 
to both the UHECRs seen by PAO, and the AGN in the GIO catalog. The 
directions are shown in an equal-area Hammer- Aitoff projection in Galactic 
coordinates; the Galactic plane is the equator (Galactic latitude b = 0°), 
and the vertically-oriented grid lines are meridians of constant Galactic lon- 
gitude, I. The star indicates the south celestial pole (SCP), the direction 
directly above Earth's south pole (effects like precession and nutation of 
the Earth's axis are negligible for this application and we ignore them in 
this description). The thick gray line bounds the PAO field of view. The 
UHECR and AGN directions are displayed as "tissots," projections of cir- 
cular patches centered on the reported directions. The small tissots show 
the UHECR directions; the tissot size is 2°, corresponding to ~ 2 standard 
deviation errors, and the tissot color indicates energy. The large green tissots 
indicate AGN directions; the tissot size is 5°, corresponding to a plausible 
scale for magnetic deflection of UHE protons in the Galactic magnetic field. ^ 
The tissots are rendered with transparency; the two darker tissots near the 
Galactic north pole indicate pairs of AGN with nearly coincident directions. 
Two of the AGN tissots are outlined in solid black; these correspond to 
the two nearest AGN, Centaurus A (Cen A, also known as NGC 5128) and 
NGC 4945, neighboring AGN at distances of 4.0 and 3.9 Mpc (as reported in 
GIO). Five others are outlined in dashed black; these have distances ranging 
from 6.6 to 10.0 Mpc. The remaining 10 AGN have distances from 11.5 to 
15.0 Mpc. Four of the AGN are outside the PAO field of view, but depend- 
ing on the scale of magnetic deflection, they could be sources of observable 
cosmic rays. 

Figure 1 shows all 69 UHECR directions; the panels in figure 2 show 
subsets of the UHECRs corresponding to the three PAO analysis periods 
described above. The thin curves (teal) show geodesies connecting each 
UHECR to its nearest AGN. In the full-catalog map, there is a notica- 



If UHECRs are comprised of heavier, more positively charged nuclei, they could suffer 
much larger deflections; see § 3.4. 
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CR Energy, 55 - 150 EeV 



Fig 1. Sky map showing directions to 69 UHECRs detected by PAO, and to 17 nearby 
AGN from the catalog of Goulding et al.. Directions are shown in an equal-area Hammer- 
Aitoff projection in Galactic coordinates. Thick gray line indicates the boundary of the 
PAO field of view. Small tissots show UHECR directions; tissot radius is 2° corresponding 
to ~ 2 standard deviation errors; tissot color indicates energy. Large green tissots indicate 
AGN directions; tissot radius is 5°. Thin curves are geodesies connecting each UHECR to 
its nearest AGN. Maps here and in Fig. 2 show UHEGR directions from different PAO 
observing periods, as labeled. 
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Fig 2. Sky maps showing directions to UHECRs detected by PAO in each of three observing 
periods (as labeled). Symbols are as in Fig. 1. 
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ble concentration of cosmic ray directions near the directions of Cen A and 
NGC 4945; a few other AGN also have conspicuously close cosmic rays. A 
concentration in the vicinity of these two closest AGN is also evident in the 
maps for periods 1 and 2 (note that the slightly extended red tissot near 
the direction of NGC 4945 corresponds to two nearly coincident UHECRs). 
Curiously, except for a single UHECR about 6° from NGC 4945, no such con- 
centration is evident in the map for Period 3, despite it having about three 
times the number of UHECRs found in earlier periods. This is a presage 
of results from our quantitative analysis that suggest the data may not be 
consistent with simple models for the cosmic ray directions, with or without 
AGN associations. 

2.4. PAO exposure. PAO is not equally sensitive to cosmic rays com- 
ing from all directions. Quantitative assessment of evidence for associations 
or other anisotropy must account for the observatory's direction-dependent 
exposure. 

Let F be the cosmic ray flux at Earth from a source at a given direction, uj, 
i.e., the expected number of cosmic rays per unit time per unit area normal 
to oj. Then the expected number of rays detected in a short time interval 
dt is FA±{t,uj)dt, where A±{t,uj) is the projected area of the observatory 
toward co at time t. The total expected number of cosmic rays is given by 
integrating over t; it can be written as Fe{uj), with the exposure map e(a;) 
defined by 



the integral is over the time intervals when the observatory was operating, 
denoted collectively by T. Supplementary Appendix A describes calculation 
of e{uj); the thick gray curves shown in the sky maps mark the boundary of 
the region of nonzero exposure. 

3. Modeling the cosmic ray data. The basic statistical problem is 
to quantify evidence for associating some number (possibly zero) of cosmic 
rays with each member of a candidate source population. The key observable 
is the cosmic ray direction; a set of rays with directions near a putative host 
comprises a multiplet potentially associated with that host. This gives the 
problem the flavor of model-based clustering (of points on the celestial sphere 
rather than in a Euclidean space), but with some novel features: 



(2) 




• The model must account for measurement error in cosmic ray proper- 
ties. 
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• Observatories provide an incomplete and distorted sample of cosmic 
rays, so the model must account for random truncation and nonuni- 
form thinning. 

• The most realistic astrophysical models imply a joint distribution for 
the properties of the cosmic rays assigned to a particular source that 
is exchangeable rather than a product of independent distributions (as 
is the case in standard clustering). 

• The number of cosmic rays is informative about the intensity scale of 
the cosmic ray sources so the binomial point process model underlying 
standard generative clustering approaches is not appropriate. 

To account for these and other complexities, we model the data using a 
hierarchical Bayesian framework with four levels: 

1. Source properties: At the top level we specify the properties of the 
sources of cosmic rays. This may include the choice of a candidate 
source population of identified objects (e.g., a particular galaxy popu- 
lation), and/or specification of the properties of a population of uniden- 
tified sources. For a given candidate source population, we must spec- 
ify source directions and cosmic ray intensities. The simplest case is a 
standard candle model, with each source having the same cosmic ray 
intensity. More generally, we may specify a (non-degenerate) distribu- 
tion of source intensities; this corresponds to specifying a "luminosity 
function" in other astronomical contexts. For a population of uniden- 
tified sources, we must specify a directional distribution (isotropic in 
the simplest case) as well as an intensity distribution. 

2. Cosmic ray production: We model the production of cosmic rays from 
each source with a marked Poisson point process model for latent 
cosmic ray properties. The incident cosmic ray arrival times have a 
homogeneous intensity measure in time, and the marks include the 
cosmic ray energies, and latent categorical labels identifying the source 
of each ray. 

3. Cosmic ray propagation: Next we model magnetic deflection of the 
rays, scattering their directions from the source directions. This re- 
quires introducing a latent arrival direction parameter for each ray. 
Here we adopt a simple phenomenological model with a single pa- 
rameter specifying a typical scattering scale between the source and 
arrival directions. As the data become more abundant and detailed, 
the framework can accommodate more complex models, e.g., with pa- 
rameters explicitly describing cosmic magnetic fields and cosmic ray 
composition. 
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Model Levels & Random Variables 

Parameters — Latent variables — Observables 



Marked Poisson 
point process for 
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Fig 3. Schematic depiction of the levels in our cosmic ray association models, identifying 
random variables appearing in each level, including parameters of interest (hold red labels), 
latent variables representing cosmic ray properties that are not directly observable (slant 
type labels), and observables (bold blue labels). 



4. Detection and measurement: Last, we model detection and measure- 
ment, accounting for truncation and thinning of the incident cosmic 
ray flux, and measurement errors for directions and energies. 

Figure 3 schematically depicts the structure of our framework, including 
identification of the various random variables appearing in the calculations 
described below. The variables will be defined as they appear in the detailed 
development below; the figure serves as visual reference to the notation. The 
figure is not a graphical model per se. Rather, our models specify probability 
distributions over a space of graphs, each graph corresponding to a possible 
set of associations of the cosmic rays with particular sources. This framework 
builds directly on an earlier multilevel Bayesian model we developed to assess 
evidence that some sources of gamma-ray bursts repeat [29]; this model, too, 
worked in terms of probability distributions over candidate assignments. See 
[28] for a broad discussion of Bayesian methods for assessing spatio-temporal 
coincidences in astronomical data. 

Our framework is designed to enable investigators to: (1) Ascertain which 
cosmic rays (if any) may be associated with specific sources with high 
probability; (2) Estimate luminosity function parameters for populations 
of astrophysical sources; (3) Estimate the proportion of all detected cos- 
mic rays generated by each population; (4) Estimate parameters describing 
the composition-dependent effects of cosmic magnetic fields; (5) Investigate 
whether cosmic rays from a single source are deflected independently or 
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share part of their deflection history (resulting in correlated deflections). 
Task (5) is not attempted here but will be investigated in the future. 

3.1. Cosmic ray source properties. We do not anticipate the UHECR 
flux passing through a volume element at the Earth to vary in time over 
accessible time scales, so we model the arrival rate into a small volume of 
space from any particular direction as a homogeneous Poisson point process 
in time. Let denote the UHECR flux from source k. is the expected 
number of UHECRs per unit time from source k that would enter a fully 
exposed spherical detector of unit cross-sectional area. A cosmic ray source 
model must specify the directions and fluxes of candidate sources. In our 
framework, a candidate source catalog specifies source directions for a fixed 
number of potential sources, Na {Na = 17 for the GIO AGN catalog). In ad- 
dition, we presume some cosmic rays may come from uncatalogued sources, 
so we introduce a background component, labeled by /c = 0, considered to be 
a population of isotropically distributed "background" sources. We presume 
the background sources to be numerous and to each have relatively low cos- 
mic ray fluxes, so that at most a single cosmic ray should be detected from 
any given background source (i.e., we do not consider clustering of cosmic 
rays assigned to the background). In this limit, the background component 
may be described by a single parameter, i<o, denoting the total flux from 
the entire background population. 

A model must specify a distribution for {.Ffc} = {i*bi F}; in astronomical 
lingo, this corresponds to specifying a "luminosity function" for the back- 
ground and source populations. As a simple starting point, we treat Fq as 
a free parameter, and adopt a "standard candle" model specifying the Na 
candidate host fluxes, F, via a single parameter as follows. We assume all 
sources emit isotropically with the same intensity, / (number of cosmic rays 
per unit time), so the flux from a source (i.e., Fk for A; > 0) can be writ- 
ten as Ffc = I /Df, (the inverse-square law), with the (known) distance 
to source k (there could also be distance- and energy-dependent attenua- 
tion due to cosmic ray-photon interactions, but the sources we consider here 
are close enough that such attenuation should be negligible). The total flux 
from the sources is Fa = Yl,k>G ^k, and we adopt Fa as the source intensity 
parameter rather than /. Thus F^ = WkFA, with the weights Wk given by 



(3) 



Wk = 




for fc = 1 to Na- 
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3.2. Top-Level Prior Specification. We must specify a prior distribution 
for Fq and F^. Earlier observations constrained the total UHECR flux. In our 
association model, the total flux is Ft = Fq + Fa- For the null model, there 
is only one top-level parameter, the total flux from an isotropic distribution 
of source directions. So we adopt Ft as a top-level parameter, common to all 
models. For association models, this motivates an alternative parameteriza- 
tion that switches from {Fq,Fa) to {Ft, /), where / = Fa/{Fq + Fa) is the 
fraction of the total flux attributed to the candidate host population. In this 
parameterization, we can specify a common total flux prior for all models. 
This is astrophysically sensible since we have results from prior experiments 
to set a scale for the total flux. It is also statistically desirable; Bayes factors 
tend to be robust to specification of priors for parameters common to models 
being compared. 

We adopt independent priors for the total flux and the associated fraction. 
If their prior densities are g{FT) and h{f), then the implied joint prior 
density for {Fq,Fa) is 



(4) 7r(Fo,FA) = — 



Fo + Fa 

where the denominator is from the Jacobian of the transformation between 
parameterizations. In general, an independent prior for Ft and / corre- 
sponds to a dependent prior for Fq and Fa- 

For the calculations below, we adopt an exponential prior with scale s for 
Ft, and a beta prior for / with shape parameters (a, 6), so 

(5) g{FT) = -e~^-/' and /,(/) = ^—/-^(l - 

where B{a,b) is the beta function. We set the hyperparameters {s,a,b) as 
follows. 

We take s = 0.01 x 4tt km~^ y~^ for all models. This scale is compatible 
with flux estimates from AG AS A and HiRes. The likelihood functions for 
Ft from those experiments are formally different from exponentials (they 
are more concentrated away from zero) , but since this prior is common to all 
models, and since the PAO data are very informative about the total flux, 
our results are very robust to its detailed specification. 

For the beta prior for /, our default choice is a = 6 = 1, which corresponds 
to a uniform prior on [0, 1] . We also repeat some computations using 5 = 5 
to investigate the sensitivity of Bayes factors to this prior. This case skews 
the prior downward, increasing the probability that / is close to 0. 
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3.3. Cosmic ray mark distributions. Given the fluxes, we model cosmic 
ray arrival times with a superposition of homogeneous Poisson point pro- 
cesses from each component. Besides its arrival time, each cosmic ray has 
a label associated with it, identifying its source component. Let A be an 
integer-valued latent label for a UHECR, specifying its source (A = for 
the background, or > 1 for AGN k). Since a superposition of Poisson 
processes is a Poisson process, we may consider the arrival times for the 
UHECRs arriving at Earth to come from a total event rate process, and the 
labels to come from a categorical mark distribution with probability mass 
function 



In the absence of magnetic deflection, the labels could be replaced by source 
directions (with background source directions assigned isotropically) , and 
the process could be considered to be Poisson in time with a directional 
mark distribution. But magnetic deflection requires a more complex setup. 

Our full framework also assigns energies as marks for each cosmic ray, 
drawn from a distribution describing the cosmic ray spectrum. The energies 
would be important in an analysis that seeks to infer the cutoff energy 
distinguishing local UHECRs rays from lower-energy cosmic rays (i.e., the 
GZK cutoff). Although the PAO-10 catalog includes event energies, the PAG 
team has already made an energy cut, and in the absence of lower-energy 
data, we cannot usefully infer a cutoff. Thus in the analysis presented here, 
we ignore the energy mark distribution. (For models with more complex 
luminosity functions and more distant sources, the energies would play a 
role in accounting for the suppression of flux from distant sources due to 
interactions with cosmic backgrounds.) 

3.4. Cosmic ray deflection. After leaving a source, UHECRs will have 
their paths deflected as they traverse galactic and intergalactic magnetic 
fields. The Galactic field is partially measured and is known to have both 
a turbulent component (varying over length scales below ~ 1 kpc) and a 
regular component (coherent over kpc scales and largely associated with 
spiral arms), with typical field strengths ~ 1 fiG. The magnetic fields of 
other galaxies are at best crudely measured and believed to be similar to 
the Galactic field. The much smaller fields in intergalactic space are only 
weakly constrained (in fact, cosmic rays might provide useful additional 
constraints); the typical field strength is probably not larger than ~ 10~^ G 
except within galaxy clusters. 



(6) 



P{\ = k\Fo,F 



Fk 
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A number of investigators have modeled cosmic ray propagation in the 
Galaxy, or in intergalactic space, using physical models based on existing 
field measurements (recent examples include [21, 20, 16, 30, 5, 24] ; see 
[34] for an overview). Roughly speaking, there are two regimes of deflection 
behavior, described here in the small-deflection limit [20]. As a cosmic ray 
with energy E and atomic number Z traverses a distance L spanning a 
regular magnetic (vector) field B, it is deflected by an angle 



ds B 



li 3 kpc 2 /iG 

where s (a vector) is an element of displacement along the trajectory; the 
field and length scales are typical for the Galaxy. If instead it traverses a 
region with a turbulent structure, with the field coherence length i <^ L, 
then the defiection will be stochastic; its probability distribution has zero 
mean, and root-mean-square (RMS) angular scale 



E \ f Srms \ f L 



-1 / R \ / r \ 1/2 / ^ X 1/2 



2.3° Z 



50EeVy V4^GyV3kpcy V50pc^ 

50EeVy vinGyVioMpcy ViMpcy 



where i?rms is the RMS field strength along the path, and quantities are 
scaled to typical galactic and intergalactic scales on the first and second 
lines, respectively. 

For a detected cosmic ray, the energy is measured fairly accurately, but 
other quantities appearing in the deflection formulae may be largely un- 
known. As noted above, there is significant uncertainty in the magnitudes 
of cosmic magnetic fields, particularly for turbulent structures. Turbulent 
length scales are poorly known. Finally, the composition (distribution of 
atomic numbers) of UHECRs is not known. Low energy cosmic rays are 
known to be mainly protons and light nuclei, but the proportion of heavy 
nuclei (with Z up to 26, corresponding to iron nuclei, the most massive sta- 
ble nuclei) increases with energy up to about 10^^ eV. At higher energies, 
inferring the cosmic ray composition is very challenging, requiring both de- 
tailed measurement of air shower properties, and theoretical modeling of the 
Z dependence of hadronic interactions at energies far beyond those probed 
by accelerators. Measurements and modeling from HiRes indicate light nu- 
clei are predominant again at ~ 1 EeV and remain so at least to ~ 40 EeV 
[35]. In contrast, recent PAO measurements indicate a transition from light 
to heavy nuclei over the range ~ 3-30 EeV [2, 11]. (The discrepancy is not 
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yet explained.) For heavy nuclei, the deflection scales in both the regular 
and turbulent deflection regimes can be large, ~ 1 rad. Some investigators 
have suggested that many or most UHECRs may be heavy nuclei originating 
from the nearest AGN, Gen A, so strongly deflected that they come from 
directions across the whole southern sky (e.g., [7, 17, 8]). 

In light of these uncertainties and the relative sparsity of UHEGRs, we 
use simple phenomenological models for magnetic deflection. In the simplest 
"buckshot" model, each cosmic ray from a particular source experiences a 
deflection that is conditionally independent of the deflection of other rays 
from that source, given a parameter, k, describing the distribution of de- 
flections. We have also devised a more complex "radiant" model that allows 
cosmic rays assigned to the same source to have correlated deflections, with 
the correlation representing a partially shared deflection history. For the 
analyses reported here, we use the buckshot model; we describe the radiant 
model further in § 6. 

The buckshot deflection model adopts a Fisher distribution for the de- 
flection angles. The model has a single parameter, k, the concentration pa- 
rameter for the Fisher distribution. The probability density for observing a 
cosmic ray from direction uj if it is assigned to source k with direction Wk is 
then 



With this deflection distribution, when a cosmic ray is generated from an 
isotropic background population, its deflected direction still has an isotropic 
distribution. Accordingly, 



The At parameter is convenient for computation, but an angular scale 
is more convenient for interpretation. The contour of the Fisher density 
bounding a region containing probability P is azimuthally symmetric with 
angular radius 9 satisfying 



(9) 



K 



exp(KW • TOfc). 



47r sinh(K) 



(10) 




(11) 




where Q denotes the cone of solid angle subtended by the contour. In plots 
showing K-dependent results, we frequently provide an angular scale axis. 
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using (11) with P = 0.683, in analogy to the "Icj" region of a normal distri- 
bution.^ 

Note that, astrophysically, k has a nontrivial interpretation. If all UHE- 
CRs are the same nuclear species (e.g., all protons), then k depends solely 
on the magnetic field history experienced by cosmic rays as they propagate 
to Earth. If UHECRs are of unknown or mixed chemical composition, then 
K conflates magnetic field history and composition. In a more complicated 
model, there could be a distribution for the values of k assigned to UHE- 
CRs (accounting for different compositions and magnetic field histories); the 
distribution could depend on source direction (accounting for known mag- 
netic field structure in the Galaxy and perhaps in intergalactic space) and 
on source distance (related to the path length in intergalactic space). 

When estimating k or marginalizing over it, we adopt a log-fiat prior 
density for k G [1,1000], 

(12) pU) = -, for 1< K < 1000. 

^ ^ ^ log 1000 K - - 

The lower limit corresponds to large angular deflection scales ~ 1 rad, such 
as might be experienced by iron nuclei. The upper limit corresponds to small 
angular deflection scales ~ 1°, such as might be experienced by protons with 
E ~ 100 EeV. 

3.5. Cosmic ray detection and measurement. Even though the arrival 
rate of UHECRs into a unit volume is constant in time in our model, the 
expected number per unit time detected from a given direction will vary as 
the rotation of the Earth changes the observatory's projected area toward 
that direction, as noted above. As a result, the Poisson intensity function 
for detectable cosmic rays varies in time for each source. 

Recall that the likelihood function for an inhomogeneous Poisson point 
process in time with rate (intensity function) r(t) has the form 

(13) eM-Nc.p)Ylr{ti)6t, 

i 

where the events are detected at times ti in detection intervals of size 5t, and 
A'exp is the total expected number in the observing interval (the integral of 
the rate over the entire observing interval). The likelihood function for the 



'in the K S> 1 limit, the Fisher density becomes an uncorrelated bivariate normal with 
respect to locally cartesian arc length coordinates about the mode on the unit sphere. The 
standard deviation of this normal distribution, 6, satisfies equation (11) with P ~ 0.683; 
for K » 1 this implies 6'^ ^ 2.30/k, or 9 ^ 86.9° /k^^^. 
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cosmic ray data has a similar form, but with adjustments due to the mark 
distribution and measurement errors. 

If the label and arrival direction for detected cosmic ray i were known, 
the factor in the likelihood function associated with that cosmic ray would 
be FkA±{uJi,ti)5t, where k = Aj. In reality, both the label and the arrival 
direction are uncertain; the PAO analysis pipeline produces a likelihood 
function for the direction to the cosmic ray, ii{uji); see equation (1). 

Introducing the uncertain direction as a nuisance parameter, with a prior 
denoted by pk{uJi\K), the likelihood factor for cosmic ray i when assigned to 
source k may be calculated by marginalizing; it may be written as Fkfk,iSt, 
with 



The cosmic ray direction measurement uncertainty is relatively small (~ 1°) 
compared to the scale over which the area varies, so we can approximate 
fk,i as 



where 9i denotes the zenith angle of UHECR i (reported by PAO-10) and 
Ai = A{ti) is the area of the observatory at the arrival time of UHECR i. 
The integral can be computed analytically; 



The total event rate for cosmic rays with the properties (direction, energy, 
and arrival time) of detected ray i combines the contributions from each 
potential source, i.e., r(tj) = Y^k^kfuA^)- 

To calculate A^exp; we must account for the observatory's exposure map. 
The effective exposure given to cosmic rays from source k throughout the 
time of the survey depends, not just on the direction to the source, but also 
on the deflection distribution, pk (and thus on k), since rays from that source 
will not arrive precisely from the source direction. The exposure factor for 
source k is 



Note that has units of area x time, and for the isotropic background 
component (A; = 0), eo{K) is a constant equal to the sky-averaged exposure 



(14) 




(15) 





(17) 
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(in the notation of Supplementary Appendix A, eo = aT/4vr). To find the 
total expected number of detected cosmic rays we sum over sources: A'exp = 

The prior probability mass function for the label of a detected cosmic ray 
is not given by (6); the terms must be weighted according to the source 
exposures. The result is 

(18) PiX. = m,F,.) = -^p^ 

We now have the ingredients needed to evaluate equation (13), gener- 
alized to include the cosmic ray marks (directions and labels) and their 
uncertainties. The resulting likelihood function is 

(19) C{Fo, F, k) = exp ^ Ffce^^ fk,,Fj}j . 

The product-of-sums factor resembles the likelihood for a finite mixture 
model (FMM), if we identify the fk^i factors as the component densities, 
and the F^ factors as the mixing weights. A common technique for com- 
puting with mixture models is to rewrite the likelihood function as a sum- 
of-products by introducing latent label parameters identifying which com- 
ponent each datum may be assigned to. In Supplementary Appendix B we 
exploit this analogy to show that the likelihood function can be rewritten 
as a sum over latent assignments of cosmic rays to sources, 

(20) £(Fo,F,.) = y: (uFr^'^^-^'A n/A.. 

A \ fc / i 

where A = {Aj} and denotes an A'^c'-dimensional sum over all possible 
assignments of cosmic rays to sources, and the multiplicity mfc(A) is the 
number of UHECRs assigned to source k according to A. We suppress the k 
dependence of ek{n) and fk,i{K) here and elsewhere to simplify expressions. 
Note that the F^ dependence (for a given A) is of the same form as a gamma 
distribution. 

Rewriting the previous expression with (Ft,/) in place of {Fq,Fa), and 
using Ffc = w^Fa = fw^Fx (for k > 1), we can rewrite C{Fq, F, k) as 

(21) C{f,FT,K) = f)"'oiX)fNc-moiX)pNc 

X 

xg-FT[(i-/).o+/E.>i-.-.] -Q^-.W-Q^^^^,. 

k>l i 
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For computations it will be helpful to have the likelihood function condi- 
tional on the label assignments, 



(22) 



P{D\X, Fo, F, k) = exp ^- ^ F^efc^ 



Nc 



n 



fXi,i 



where k runs over the host labels (from to N^), and i runs over the UHECR 
labels (from 1 to Nq)- We can recover the likelihood for Fq, F, and k, by 
multiplying by the prior for A from equation (18) and marginalizing, giving 
equation (20). 

3.6. Estimating k. To estimate the deflection parameter, k, we need the 
marginal likelihood Cm{K) = P{D\k) = f dFr f dfP{D,FTj\K). The inte- 
grand is the product of equation (21) and the flux priors. Using the exponen- 
tial and beta priors described above, we have that the marginal likelihood 
for K is 



■^m(^) 



E 



(23) 



sB{a,b) 

j!Nc-mo{X)+a~l^-^ 



f) 



mo(A)+6-l 



i + (l-/)6o + /Efc>lW^fcefc 



Nc+l 



df. 



Computing Cmin) requires summing over all possible values of A which is 
intractable in practice. In Supplementary Appendix C.2, we describe how 
to use Chib's method [12] to calculate this marginal likelihood. 

3.7. Model Comparison. We use Bayes factors (ratios of marginal likeli- 
hoods) to compare models both conditioned on k (using marginal likelihood 
functions Cm{i^)), and after marginalizing over k (using the log- flat prior of 
equation (12), and numerical quadrature over k). 

We consider three models. The null model, Mq, assumes that all the UHE- 
CRs come from the isotropic background source population; recall that it has 
no K dependence (see equation (10)). Model M\ allows the UHECRs to come 
from any of the 17 AGN in the catalog or from the isotropic background. We 
also consider another model, M2, in which the UHECRs may come from the 
isotropic background or either of the two closest AGN, Cen A (NGC 5128) 
and NGC 4945; this model is motivated in part by recent suggestions that 
most UHECRs may be heavy nuclei from a single nearby source, as cited 
above. (We also briefly explore a similarly-motivated fourth model that as- 
signs all UHECRs to Cen A; as noted below, this model is tenable only for 
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K ~ 0.) In order to compare models Mi and M2 (conditioned on k) to the 
null model, we compute the Bayes factors: 



where Cm,o is the marginal likelihood for the null model (similar equations 
hold for models that marginalize over k). The value of Cmfi can be found 
from equation (20), noting that for the null model, there is only one term in 
the sum over A (with all Aj = 0, since the only allowed value of /c is /c = 0). 
Marginalizing this term over Ft (equal to in this case) gives 



3.8. Computational techniques. The principal obstacle to computing with 
this framework is the combinatorial explosion in the number of possible as- 
sociations as the sizes of the candidate source population and the cosmic ray 
sample grow. For small amounts of magnetic deflection, the vast majority 
of candidate associations are improbable (they associate well-separated ob- 
jects with each other). But there is evidence that UHECRs may be massive 
(and thus highly charged) nuclei, which would undergo significant deflec- 
tion. To probe the full variety of astrophysically interesting models requires 
techniques that can handle both the small- and large-deflection regimes, for 
catalog sizes corresponding to current and forthcoming catalogs from PAO. 

For parameter estimation within a particular model, we have developed 
a Markov chain Monte Carlo (MCMC) algorithm that draws samples of 
the parameters /, Ft, and A from their joint posterior distribution. The 
algorithm takes advantage of two features of the models described above. 
First, by introducing latent labels, A, we could write the likelihood function 
in a sum-of-products form, equation (20), with factors that depend on the 
fluxes F^ in the manner of a gamma distribution. Second, the forms of the 
likelihood and priors are conjugate for Ft and the labels, so we can find 
closed-form expressions for their conditional distributions. These features 
enable us to use Gibbs sampling techniques well known in mixture modeling 
for sampling the Ft and A parameters. We handle the / parameter using a 
random walk Metropolis algorithm, so our overall algorithm is a Metropolis- 
within-Gibbs algorithm. Supplementary Appendix C provides details on its 
implementation. 

We treat the deflection parameter, k, specially, considering a logarithmically- 
spaced grid of values that we condition on. We did this so that we could 
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explore the k dependence more thoroughly than would be possible with 
posterior sampling of k. Of course, our Metropolis-within-Gibbs algorithm 
could be supplemented with k proposals to enable sampling of the full pos- 
terior. 

Finally, using Bayes factors to compare rival models requires computing 
marginal likelihoods, which are not direct outputs of MCMC algorithms. Us- 
ing a simplified version of our model, and modest-sized simulated datasets, 
we explored several approaches for marginal likelihood computation in a 
regime where we could compute the correct result via direct summation over 
all feasible associations. We explored the harmonic mean estimator (HME), 
Chib's method, and importance sampling algorithms. The HME performed 
poorly, often apparently converging to an incorrect result (such behavior 
is not unexpected; see [39]). Importance sampling proved inefficient. Chib's 
method was both accurate and efficient in these trial calculations, and be- 
came our choice for the final implementation. Supplementary Appendix C 
provides details. 

4. Results. Recall that the UHECR data reported by PAO-10 are di- 
vided into three periods. The PAO team used an initially larger period 1 
sample (including lower-energy events) to optimize an energy threshold de- 
termining which events to analyze in period 2; the reported period 1 events 
are only those with energies above the optimized threshold. The optimiza- 
tion maximized a measure of anisotropy in the above-threshold Period 1 
sample. Without access to the full Period 1 sample, we cannot evaluate the 
impact of this optimization on our modeling of anisotropy in the reported 
Period 1 data (nor can we usefully pursue a Bayesian treatment of a GZK 
energy cutoff parameter). Because of this complication, we have performed 
analyses for various subsets of the data. As our main results, we report cal- 
culations using data from Periods 2 and 3 combined ("untuned data") and 
for Periods 1, 2, and 3 combined ("all data"). We also report some results 
for each period considered separately, and we use them to perform a simple 
test of consistency of the results across periods, in an effort to assess the 
impact of tuning on the suitability of the Period 1 data for straightforward 
statistical analysis. 

4.1. Results conditioning on the deflection scale, k. We first consider 
models conditional on the value of the magnetic deflection scale parame- 
ter, K, calculating Bayes factors comparing models, and estimates of the 
association fraction, /. 

We report model comparison results as curves showing Bayes factors as 
functions of k. These quantities are astrophysically interesting but must be 
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interpreted with caution. The actual values of the conditional Bayes factors 
can only be interpreted as Bayes factors for a particular value of k deemed 
interesting a priori. For example, were one to assume that UHECRs are 
protons, adopt a particular Galactic magnetic field model, and assume that 
intergalactic magnetic fields do not produce significant deflection (which is 
plausible for protons from local sources), one would be interested only in 
large values of k of order several hundred (corresponding to small angular 
scales for deflection). On the other hand, if one presumed that UHECRs 
are predominantly heavy nuclei, then deflection by Galactic fields could be 
very strong, corresponding to k of order unity (deflection by intergalactic 
fields might also be significant in this case) . Models hypothesizing that most 
UHECRs are heavy nuclei produced by Cen A would fall in this small-K 
regime. By presenting results conditional on k, various cases such as these 
may be considered. Also, the Bayes factor conditioned on n is proportional 
to the marginal likelihood for k, so the same curves summarize the infor- 
mation in the data for estimating n if it is considered unknown. We plot 
the curves against a logarithmic k axis, so they may be interpreted (up to 
normalization) as posterior probability density functions based on a log-flat 
K prior. 

The Bayes factors comparing models Mi and M2 to Mq for various values 
of At G [1, 1000], and for various partitions of the data, are shown in Figure 4. 
For cases using only the untuned data (periods 2, 3, and 2-1-3), we find 
that both BFio and BF20 (see equation 24) are close to 1 for all values 
of At € [1, 1000] for the beta-(l,l) (uniform) prior for /. The Bayes factors 
are only a little higher in the case of beta-(l,5) prior, indicating the results 
are robust to reasonable changes in the / prior. These values imply that the 
posterior odds for the association models Mi and M2 versus the null isotropic 
background model Mq are nearly equal to the prior odds, indicating the 
untuned data provide little evidence for or against either association model 
versus the isotropic model. 

Considering the Period 1 data qualitatively changes the results. The solid 
(blue) curves in Figure 4 show the Bayes factor vs. k results based solely on 
the Period 1 data; there is strong evidence for association models conditioned 
on K values of around 50 to 100.^*^ Analyzing the data from all three periods 
jointly produces the long-dashed (purple) curves. Using a uniform prior for 
/, we find BFio attains a maximum of 90 at k ~ 46 while BF20 attains a 
maximum of 262 at k ~ 38. Both BFio and BF20 are larger than 30 for all 

^"A common convention for interpreting Bayes factors is due to Kass and Raftery, who 
consider a Bayes factor between 3 and 20 to indicate "positive" evidence and between 20 
and 150 to indicate "strong" evidence [25]. 
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Fig 4. Bayes factors comparing the assocation model with 17 AGN (top row) or 2 AGN 
(bottom row) with the null isotropic background model, conditional on k, shown as a func- 
tion of K (bottom axis) and the corresponding deflection angle scale, a (top axis). Results 
are shown for various partions of the data (identified by line style, identified in the legend), 
and for two choices of the prior on f : a flat prior (left column), and a Beta{l, 5) prior 
(right column). 
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K G [20, 120]. Both of the association models are strongly preferred over the 
null in this range of k, while the comparison is inconclusive for k outside 
this range. 

The originally published data (in PAO-08) covered periods 1 and 2. For 
comparison with studies of that original catalog, Figure 4 include curves 
showing the Bayes factor vs. k based on data from periods 1 and 2. This 
partition of the data produces the largest Bayes factors, ~ 1000 for k k. 50. 
The curves are qualitatively consistent with accumulation of evidence from 
Period 1 and Period 2}^ These results amplify what was found in the analysis 
using all of the data: the strongest evidence for association comes from the 
Period 1 data. This is troubling because this data was used (along with 
unreported lower-energy data) to tune the energy cut defining all of the 
samples, and there is no way for independent investigators to account for 
the effects of the tuning on the strength of the evidence in the Period 1 data. 

We show marginal posterior distributions for / in Figure 5, for both M\ 
and M2, using both the untuned data, and using all data. For a given model, 
the posterior does not change much when Period 1 data are included. The 
posteriors indicate evidence for small but nonzero values of /, of order a few 
percent to 20%. They strongly rule out values of / > 0.3, indicating that 
most UHECRs must be assigned to the isotropic background component in 
these models. This holds even for values of k as small as ~ 10, corresponding 
to quite large magnetic deflection scales, as might be experienced by iron 
nuclei in typical cosmic magnetic fields. Of course, when k = the associa- 
tion models become indistinguishable from an isotropic background model. 
These results suggest that a model assigning all UHECRs to a single nearby 
source, such as Cen A, may be tenable only with very large deflection scales. 
In Supplementary Appendix D we briefly explore such a model and confirm 
this conclusion. 

A recent approximate Bayesian analysis [38] , based on a discrete pixeliza- 
tion of the sky, attributed a similar fraction of the sample of 27 Period 1 
and Period 2 UHECRs to standard-candle AGN sources, considering ~ 900 
AGN within 100 Mpc from the VCV as candidate sources. We compare our 
approaches and results in Supplementary Appendix E. 

The posterior mode is at larger values of / for model Mi (with 17 AGN) 
than for M2 (with the two closest AGN), suggesting that there is evidence 
that AGN in the GIO catalog besides Cen A and NGC 4945 are sources of 

^^Note that the Bayes factor for the 1 + 2 partition should not be expected to equal 
the product of the Bayes factors based on the Period 1 and Period 2 partitions, because 
the models are composite hypotheses, and the data from different periods generally will 
favor different values of the model parameters. 
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UHECRs. Our multilevel model allows us to address source identification 
explicitly, by providing a posterior distribution for possible association as- 
signments (values of A) . In Table 1 we show marginal posterior probabilities 
for associations that have non-negligible probabilities (i.e., > 0.1), based on 
models Mi and M2 for two representative values of k (k = 31.62, corre- 
sponding to a 15.5° deflection scale, is a favored value for analyses including 
Period 1 data as shown below; k = 1000, corresponding to a 2.7° deflection 
scale, may be appropriate if UHECRs are predominantly protons). Rows 
are labeled by cosmic ray number, i, and columns by AGN number, k; the 
tabulated values are P{Xi = k\- ■ ■). Cosmic rays 17 and 20 (in Period 2) 
are associated with Cen A (AGN 13) with modest to high probability in 
all cases. No other assignments are robust (notably, Period 3 has no robust 
assignments, despite containing more than three times the number of cos- 
mic rays as Period 2). If UHECRs experience only small deflections, then 
besides the two Cen A associations, it is highly probable that cosmic ray 8 
(in Period 1) is associated with NGC 4945. For the larger deflection scale, 
nearly a quarter of the cosmic rays have candidate associations with prob- 
ability > 0.1, although none of those associations have probability > 0.5. 
The larger favored value of / for Mi thus reflects the 17 AGN model find- 
ing enough plausible associations (besides those with Cen A and NGC 4945) 
that it is likely that some of them are genuine, even though it cannot specify 
which. 

We can also calculate posterior probabilities for multiplet assignments. 
In general, the probability for a multiplet assigning a set of cosmic rays to 
a particular candidate source will not be the product of the probabilities 
for assigning each ray to the source. In Table 1 we see that CRs 17 and 
20 are often commonly assigned to Cen A. As an example, for Mi with 
K = 1000, their separate probabilities for assignment to Cen A are 0.85 and 
0.94, respectively. The probability for a doublet assignment of both of them 
to Cen A in this model is 0.80, which happens to be nearly equal to the 
product of their separate (marginal) assignment probabilities. Were we to 
marginalize over k, the multiplet probability would differ from the product, 
since the preferred value of k differs slightly between these two CRs. 

4.2. Results with n as a free parameter. Joint marginal posterior dis- 
tributions for \ogiQ(K,) and / are shown in Figure 6, for both association 
models, and for untuned data and all data samples. For the all-data cases, 
the joint posterior distribution is unimodal and attains its maximum at 
(k=32 , /=0.13) and {k=32 , /=0.09) for the association model with 17 
AGN and 2 AGN, respectively. For untuned data, the joint posteriors are 
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Table 1. The posterior probability that each cosmic ray is assigned to each AGN given n — 31.62 and 1000, using cosmic rays from 
periods 1+2+3. Only assignments with probabilities greater than 0.1 are shown. The AGN identifiers are: 2: NGC 0613; 7; NGC 3621; 
11: NGC 4945; 13: NGC 5128 (Cen A); 17: NGC 6300. 
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Fig 6. Marginal joint posterior distributions for the magnetic deflection concentration 
parameter, k, and the association fraction, f, considering UHECR data from different 
periods, and candidate host catalogs of 2 or 17 nearby AGN. Contours bound HPD credible 
regions of probability 0.25 (blue), 0.5 (green), 0.75 (red), 0.95 (brown), and 0.99 (gray). 



bimodal with one of the mode at the value of k shghtly less than in the 
case of all 3 periods, and the other mode at k ~ 1000, similar to the plot 
of Bayes factors in Figure 4. The results from the two samples are more 
similar than this description may indicate; they have significant peaks in 
the same region, but the likelihood function is relatively flat for the largest 
and smallest values of k (this is also apparent in Figure 4). 

In all cases, the preferred values of k correspond to deflection scales w 10°. 
As noted above, models of proton propagation in cosmic magnetic fields pre- 
dict deflections of a few degrees. The posterior distributions for k are com- 
fortably consistent with such predictions, but they do favor the larger scales 
that would be experienced by heavier nuclei. These scales are consistent 
with the suggestive evidence from PAO that UHECRs may be comprised of 
heavier nuclei than lower-energy cosmic rays. 

Values for Bayes factors accounting for k uncertainty are listed in Table 2, 
for both association models, and for both individual and combined data 
samples (these values are based on the default flat prior for /). We find 
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Table 2 

Overall Bayes factors comparing association models with 17 AGN or 2 AGN to the null 
isotropic background model, for two different priors for f 

strong evidence for both association models when considering all the cosmic 
ray data. If we exclude the tuned data of period 1, then we see positive 
evidence for association if we consider only period 2 but positive evidence 
for the null model if we consider only period 3. If we pool the untuned 
data, the data are equivocal. Together, these results raise concerns about 
consistency of the data and adequacy of the models; we address this further 
below. These results do not change qualitatively when we use the alternative 
prior for / described in § 3.2. 

Marginal posterior distributions for / and for Ft are shown in Figure 7. 
For the untuned data, the posterior mode of / is 0.051 for Mi (17 AGN) and 
0.047 for M2 (2 AGN); the 95% highest density credible intervals for / are 
[0,0.23] and [0.002,0.145], respectively. Using all of the data, the distribu- 
tions shift to somewhat larger values of /; the posterior mode of / is 0.11 for 
Ml and 0.08 for M2, and / = has a significantly smaller density. However, 
the uncertainties are large enough that the / estimates are consistent with 
each other. The posterior distributions for Ft are very similar in all models. 
The peaks are a little higher and the widths of the peaks are smaller when 
we consider the cosmic rays from periods 1-3, as expected, since we have 
more data. The posterior modes correspond to total fluxes of about 0.04 
km~^ yr""*^ in all cases. 

5. Model checking. In this section we describe results of two types of 
tests of our models. Both are motivated by the evident variability of some 
of the parameter estimation and model comparison results presented above 
with the choice of observing period or periods to include in the analysis. 

Our first tests address whether the variability indicates that the properties 
of the detected cosmic rays change from period to period, presuming that 
one of our models can adequately describe the data within each period. 
We implement a simple Bayesian change-point analysis that shows there 
is no significant evidence for variability of model parameters from period 
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Fig 7. Marginal posterior distributions for f (the fraction of UHECRs associated with 
AGN in candidate catalogs), and Ft (the total flux), considering UHECR data from dif- 
ferent periods, and models associating UHECRs with either 2 or 17 nearby AGN. 



to period. That is, presuming one of the models is adequate, the apparent 
discrepancy among the Bayes factors in Table 2 reflects variability that may 
be expected for these modest sample sizes. 

Alternatively, we may consider the possibility that none of the considered 
models accurately describes the data, in which case the variability could be 
an indication of incompatibility of the data with all of the models. To address 
this, we perform predictive checks using the Bayes factors for subsamples of 
the data as test statistics: we ask whether the period-to-period Bayes factor 
variations we find using the observed data are surprising compared with 
what one finds from simulations under various models. We first compare 
the observed Bayes factors with predictions from the null model, and then 
from representative association models. These tests are meant to explore 
broad compatibility of predictions and observations; we make no attempt to 
formally assign significances to the comparison such as p-values. 

5.1. Change-point analysis. As seen above, including the data from Pe- 
riod 1 can change Bayes factors dramatically (but not estimates of Ft, /, or 
k); also, Bayes factors differ markedly even between the untuned samples. 
Periods 2 and 3. Based on astrophysical considerations, the properties of 
incident UHECRs should not vary over the time scales under consideration. 
Even if CRs are generated in brief bursts, these bursts would have observed 
durations of hundreds or thousands of years because the cosmic rays would 
take paths of varying lengths to Earth. Therefore, any evidence indicating 
that the observed properties of cosmic rays are changing over a period of a 
few years would indicate a problem with the data, e.g., statistical inhomo- 
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geneity due to the special treatment of data in period 1, or instability of the 
observatory's apparatus or data reduction pipeline. 

As a simple test for variability in the ensemble properties of cosmic rays 
from period to period, we compared versions of Mi and M2 that allow 
model parameters to change between periods to versions that keep the pa- 
rameters the same for all periods. We describe implementation and results 
of this change-point analysis in Supplementary Appendix F. It shows that 
the period-to-period variation of Bayes factors for association vs. isotropy is 
consistent with no variation of cosmic ray properties across period bound- 
aries. 

5.2. Null model predictive checks. If the null model is in operation, a 
measure of surprise for a particular data set would be a high Bayes factor 
favoring an association model. Of particular concern here are the large Bayes 
factor values we find for the Period 1 data, which may not be representative 
because of tuning. We ask: under the null isotropic background model and 
in the absence of tuning, how likely it is to see Bayes factors as high as 100 
for some values of k in Period 1 (as shown in Figure 4)? 

To address this, we generate 200 datasets, each of which has 14 CRs (the 
size of the Period 1 sample). The CR directions are generated uniformly 
over the sky, and are accepted with probability proportional to the expo- 
sure map for that direction. We calculated Bayes factors for the 17 AGN 
association model with k = 10,31.62,100,316.2 and 1000. Figure 8 shows 
cumulative histograms (as survival functions) of the resulting Bayes factors, 
for each value of k. Most of the datasets have Bayes factors less than 1. The 
smallest Bayes factor we found is 0.04. Out of these 200 datasets, we found 
only 3 datasets with Biq > 10. Each of these datasets has Biq > 10 only 
at one value of k that we computed. The Bayes factors are 29 (k = 10), 
44 (k = 31.62) and 13 {k = 316.2). We conclude that the Bayes factor 
greater than 100 seen in Period 1 is unlikely to be due to chance if UHECRs 
are fairly sampled from isotropically distributed directions. This implies the 
distribution of directions in the Period 1 sample is anisotropic, but the cal- 
culation does not address whether this may be due to tuning or to genuine 
anisotropy. 

5.3. Association model predictive checks. Now we address compatibility 
of the discrepant Bayes factors with the predictions of the association mod- 
els. We ask: under the association model with 17 AGN, how likely is it to 
simultaneously see large Bayes factors in periods 1 or 2 (which have small 
numbers of detected CRs) and a small Bayes factor in Period 3 (which has 
a larger number of detected CRs)? 
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Fig 8. The number of simulated datasets having Bayes factor > B vs. B, based on 200 
simulated datasets with 14 CRs generated under the null isotropic background model. Axes 
are logarithmic. 



To address this, we generate 100 datasets, each with 69 CRs, partitioned 
into samples of size 14, 13 and 42 for periods 1, 2 and 3, respectively. We 
simulate from a model with association fraction / = 0.1, near the modes 
found in our analyses of the PAO data. We first simulate an incident flux of 
cosmic rays by assigning a ray to the AGN population with probability equal 
to /, or to the isotropic background with probability 1 — /. AGN-generated 
CRs get assigned to one of the 17 AGN with probabilities proportional to 
the inverse square distances of the AGN. The arrival directions of these 
CRs are drawn from a Fisher distribution centered at the source AGN, 
with K = 50. Each generated CR direction is accepted with probability 
proportional to the exposure map for that direction; a dataset is generated 
when 69 candidate events are accepted. Subsets of each simulated dataset, 
with sizes corresponding to those of the PAO subsamples, were analyzed 
with models conditioning on various values of k. Figure 9 shows cumulative 
histograms of the resulting Bayes factors. 

Out of 100 simulated datasets, we found 17 that have some of the Bayes 
factors > 30. Of these, 9 datasets have Biq > 30 in either Period 1 or 2, 
and a low Biq value in Period 3, while the other 8 datasets have -Bio > 30 
in Period 3 and low Biq values in periods 1 and 2. Out of these 17 datasets, 
10 of them have some values of Biq > 100; two have Bayes factors over 
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Fig 9. The number of simulated datasets having Bayes factor > B vs. B, based on 100 
simulated datasets with CRs generated under the 17-AGN association model with f — 0.1 
and K = 50. Panels show results for datasets with sizes corresponding to the PAO samples 
for periods 1-3, as labeled; colors of curves indicate the k value used for analyzing the 
simulated data. Axes are logarithmic. 



1000 (in either Period 1 or Period 2). These results indicate that Bayes 
factors as high as 100 (as seen in Period 1) may be expected a few percent 
of the time under the association model, and also that large between-period 
variations in the Bayes factors should be expected. However, since Bayes 
factors >10 are uncommon, and since the distribution is independent from 
one interval to the next, it is unlikely to see significantly large Bayes factors 
in both periods 1 and 2. We crudely estimate the probability for seeing Bayes 
factors whose product is > 1000 as less than ~ 10"'^. Thus the pattern of 
Bayes factors is surprising, again indicating that the tuning of the Period 1 
data is problematic. 

6. Summary and Discussion. We have described a new multilevel 

Bayesian framework for modeling the arrival times, directions, and energies 
of UHECRs, including statistical assessment of directional coincidences with 
candidate sources. Our framework explicitly models cosmic ray emission, 
propagation (including deflection of trajectories by cosmic magnetic fields), 
and detection. This approach cleanly distinguishes astrophysical and exper- 
imental processes underlying the data. It handles uncertain parameters in 
these processes via marginalization, which accounts for uncertainties while 
allowing use of all of the data (in contrast to hypothesis testing approaches 
that optimize over parameters, requiring holding out a subset of the data for 
tuning) . We demonstrated the framework by implementing calculations with 
simple but astrophysically interesting models for the 69 UHECRs with ener- 
gies above 55 EeV detected by PAO and reported in PAO-10. Here we first 
summarize our findings based on these models, and then describe directions 
for future work. 
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6.1. Astrophysical results. We modeled UHECRs as coming from either 
nearby AGN (in a volume-limited sample including all 17 AGN within 
15 Mpc) or an isotropic background population of sources; AGN are consid- 
ered to be standard candles in our models. We thoroughly explored three 
models. In Mq all CRs come from the isotropic background; in Mi all CRs 
come from either a background or one of the 17 closest AGN; in M2 all CRs 
come from either a background source or one of the two closest AGN (Gen A 
and NGG 5128, neighboring AGN at a distance of 5 Mpc). The data were 
reported in three periods. Data from Period 1 were used to tune the energy 
threshold defining the published samples in all periods by maximizing an 
index of anisotropy in Period 1. Out of concern that this tuning compro- 
mises the data in Period 1 for our analysis, we analyzed the full dataset and 
various subsamples, including an "untuned" sample omitting Period 1 data. 

Using all of the data, Bayes factors indicate there is strong evidence fa- 
voring either Mi or M2 against Mq but do not discriminate between Mi 
and M2. The most probable models associate about 5% to 15% of UHECRs 
with nearby AGN, and strongly rule out associating more than w 25% of 
UHECRs with nearby AGN. Most of the high-probability associations in the 
17 AGN model are with the two closest AGN. 

However, if we use only the untuned data, the Bayes factors are equivocal 
(although the most probable association models resemble those found using 
all data). If we subdivide the untuned data, we find positive evidence for 
association using the Period 2 sample, but weak evidence against association 
using the much larger Period 3 sample. Together, these results suggest that 
the statistical character of the data may differ from period to period, due 
to tuning of the Period 1 data or other causes. 

One way to explore this is to ask whether the data from the various periods 
are better explained using models with differing parameter values rather 
than a shared set of values. We investigated this via a change-point analysis 
that considered the time points bounding the periods as candidate change 
points. The results are consistent with the hypothesis that the parameters 
do not vary between periods, justifying using the combined data for these 
models. This suggests the variation of the Bayes factors across periods is a 
consequence of the modest sample sizes. However, the change point analysis 
does not address the possibility that none of the models is adequate, with 
model misspecification being the cause of the apparently discrepant Bayes 
factors. 

We used simulated data from both the isotropic model and high-probability 
association models to perform predictive checks of our models, using the 
Bayes factors based on subsets of the data as test statistics. Simulations 
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based on the isotropic model indicate that large Bayes factors favoring asso- 
ciation are unlikely for untuned samples of the size of the Period 1 sample. 
Simulations based on representative association models indicate that such 
Bayes factors are not surprising for samples of the size of Period 1, consid- 
ered in isolation. But the observed pattern of large Bayes factors for the 
subsamples in periods 1 and 2, and a small Bayes factor for the much larger 
Period 3 subsample, is very surprising. The full dataset thus is not fit com- 
fortably by either isotropic models or standard-candle association models. 
Whether the effects of tuning could explain the apparent inconsistencies re- 
mains an open question that is not easy to address without access to the 
untuned data. 

Restricting to the untuned data (periods 2 and 3), the pattern of Bayes 
factors is consistent with both isotropic models and representative standard 
candle association models. The best-fitting association models assign a few 
percent of UHECRs to nearby AGN; at most ~ 20% may be associated 
with AGN, with the remainder assigned to sources drawn from an isotropic 
distribution. Magnetic deflection angular scales of ~ 3° to 30° are favored. 
Models that assign a large fraction of UHECRs to a single nearby source 
(e.g., Gen A) are ruled out unless very large deflection scales are specified a 
priori, and even then they are disfavored. 

Even restricting to results based on the untuned data, we hesitate to offer 
these models as astrophysically plausible explanations of the PAG UHEGR 
data, both because of how important the problematic Period 1 sample is in 
the analysis, and because of astrophysical limitations of the models consid- 
ered here and elsewhere. In particular, the high-probability models assign 
the vast majority of UHEGRs to sources in an isotropic distribution. But the 
observation by PAG of a GZK-like cutoff in the energy spectrum of UHE- 
GRs argues strongly that UHEGRs originate from within ~ 100 Mpc, where 
the distribution of both visible matter (galaxies) and dark matter is signif- 
icantly anisotropic. If most or all UHEGRs are protons, so that magnetic 
deflection is not very strong, an isotropic distribution of UHEGR arrival 
directions is implausible. It then may be the case that some of the strength 
of the evidence for association with nearby AGN is due to the "straw man" 
nature of the isotropic alternative. On the other hand, if most UHEGRs 
are heavy nuclei, then strong magnetic deflection could isotropize the ar- 
rival directions. The highest probability association models have relatively 
small angular deflection scales, but it could be that the few UHEGRs that 
these models associate with the nearest AGN happen to be protons or very 
light nuclei. Future models could account for this by allowing a mixture of 
K values among cosmic rays, as noted in § 3.4. 
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In addition, the standard candle cosmic ray intensity model adopted here 
and in other studies very likely strongly constrains inferences. The strongest 
visible clustering of measured UHECR directions is toward the two closest 
AGN, just 5 Mpc away and only a few degrees apart on the sky. Most of 
the high-probability associations identified in our models are to these AGN. 
They are so close that they imply standard-candle cosmic ray intensities 
that would produce a negligible flux of cosmic rays from the vast majority 
of other AGN within 100 Mpc. Put another way, a standard-candle model 
assigning just one UHECR to an AGN near the 100 Mpc GZK limit would 
imply a cosmic ray flux from nearby AGN so huge that this scenario is ruled 
out simply by visible inspection of sky maps. Models with more flexible 
luminosity functions would likely allow assignment of many more UHECRs 
to sources spanning a range of distances. 

6.2. Future directions. All of these considerations indicate a more thor- 
ough exploration of UHECR production and propagation models is needed. 
We thus consider the analyses here to be a demonstration of the utility and 
feasibility of analyzing such models within a multilevel Bayesian framework, 
and not a definitive astrophysical analysis of the data. We are pursuing 
more complex models elsewhere, expanding on the present analysis in four 
directions. 

First, we are considering larger, statistically well-characterized catalogs 
of potential hosts, e.g., the recently-compiled catalog of X-ray selected AGN 
detected by the Burst and Transient (BAT) instrument on the Swift satellite, 
a catalog considered by PAO-10. 

Second, we are building more realistic background distributions, for ex- 
ample by using the locations of nearby galaxy clusters, or the entire nearby 
galaxy distribution, to build smooth background densities (e.g., via kernel 
density estimation, or fitting of mixture or multipole models). 

Third, we are considering richer luminosity function models, including 
models assigning a distribution of cosmic ray intensities to all candidate 
sources, and models that place some sources in "on" states and the others 
"off." The latter models are motivated both by the possibility of beaming of 
cosmic rays, and by evidence for AGN intermittency in jet substructure, and 
could enable assignment of significant numbers of UHECRs to both distant 
and nearby sources. 

Finally, more complicated deflection models are possible. For example, 
we have developed a class of "radiant" models that produce correlated de- 
flections (as seen in some astrophysical simulations). For a radiant model, 
each source has a single guide direction associated with it, drawn from a 
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Fisher distribution centered at the source direction, with concentration Kg] 
the guide direction serves as a proxy for the shared magnetic deflection 
history of cosmic rays from that source. Each cosmic ray associated with 
that source then has its arrival direction drawn from an independent Fisher 
distribution centered about the guide direction, with concentration poten- 
tiaUy depending on cosmic ray energy and source distance; this distribution 
describes the effect of the deflection history unique to a particular cosmic 
ray. The resulting directions for a multiplet will cluster along a ray pointing 
toward the source. The resulting joint distribution for the directions in a 
multiplet (with the guide direction marginalized) is exchangeable but not 
independent. 

For the current, modest-sized UHECR catalog, the complexity of some 
of these generalizations is probably not warranted. But PAO is expected 
to operate for many years, and the sample is continually growing in size. 
Making the most of existing and future data will require, not only more 
realistic models, but also more complete disclosure of the data. In particular, 
a fully Bayesian treatment — including modeling of the energy dependence in 
the UHECR flux and deflection scale — requires data uncorrupted by tuning 
cuts. Further, the most accurate analysis should use event-specific direction 
and energy uncertainties (likelihood summaries), rather than the typical 
error scales currently reported. We hope our framework helps motivate more 
complete releases of future PAO data. 

APPENDIX A: AUGER OBSERVATORY EXPOSURE 

PAO can reliably detect and measure UHECRs arriving from directions 
within 60° from the observatory zenith. Due to Earth's rotation, the zenith 
traces a circular path on the sky, and the PAO field of view changes with 
time. The observatory's geodetic latitude is —35.5°, so the field of view al- 
ways includes the SCP and directions within a 5.5° cone about it. Outside 
of that cone, the time spent within the field of view decreases with increas- 
ing latitude, vanishing for northern latitudes above 24.5°. This boundary 
corresponds to the thick gray curve shown in the sky maps. 

In addition, at a given instant, the projected area of the observatory varies 
with direction. Since the observatory detects air showers, the effective area 
of the detector toward a particular source direction is the projected area of 
the layer of atmosphere above PAO toward that direction. Due to Earth's 
rotation, this projected area is a function of time; we denote it by A±{t,uj) 
for direction co (a unit vector toward a fixed direction on the sky) at time t. 
We account for the zenith angle criterion by setting A± = for directions 
outside of the instantaneous field of view. 



MULTILEVEL MODELS FOR COSMIC RAYS 



43 



0.6 




-90 -60 -30 30 60 90 

Declination (deg) 

Fig 10. Average projection factor, m{uj), describing the declination dependence of the PAO 
exposure map. 

As described in § 2.4, the exposure map, e(a;), is defined by 
(26) e{uj) := J A±{uj,t)dt. 

The projected area can be written as A±{io,t) = A{t) fi{tj , t) , where A{t) is 
the effective planar area of the detection volume, and is a projection 

factor. A{t) varies slowly with time as the observatory grows. The projec- 
tion factor varies much more quickly, due to Earth's rotation. The PAO 
team has shown that for UHECRs, to a good approximation a simple geo- 
metric projection varying periodically with a period of 1 sidereal day gives 
a very accurate description of the PAO exposure [36]. As a result, the time 
integral in equation (26) can be approximated as J A{t)m{uj)dt, where m(a;) 
is the geometric projection factor averaged over 1 sidereal day. The average 
projection factor is constant with respect to right ascension (the equatorial 
sky coordinate corresponding to geodetic longitude) due to rotational av- 
eraging. It varies strongly with declination (the equatorial sky coordinate 
corresponding to geodetic latitude). Figure 10 shows the average projection 
as a function of declination. 

With these approximations, to evaluate e{uj), we need the time integral 
of the observatory area, A{t). By convention, this quantity is reported indi- 
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rectly by describing the observatory's sensitivity to an isotropic distribution 
of sources (lower-energy cosmic rays have an isotropic distribution) . For such 
a distribution, the expected number of rays would be proportional to the 
total exposure, the integral of the exposure map over the whole sky: 

(27) ar ■= J e{uj)duj = J A{t)dt J m{oj)duj. 

The total exposure has units of area x time x solid angle; (it has also been 
called "aperture," as "exposure" is more traditionally used for quantities 
with units of area x time). The PAO team reports ax for each observing 
period in PAO-10 [33]. 

To calculate the exposure map from ut and m(a;), define the integrated 
projection factor by M := J dujm{u) (with units of solid angle). Then the 
exposure is 

(28) eiu) = ^m{co). 

APPENDIX B: LIKELIHOOD FACTORIZATION 

Consider estimating the probability density function for a sample of N 
points {xi} in a Euclidean space using a mixture model built from the pa- 
rameterized kernel density f{x; 6) with parameters 6 (typically including at 
least location and scale parameters). With the number of components in the 
mixture fixed as K, the likelihood function for the set of kernel parameters 
{9k} and (normalized) mixing weights {w^} is 



N 



(29) £{{wk},{Ok}) = ll 



i=l 



K 



^Wkfixi^Ok 



.k=l 



The factor associated with a particular datum, ^i.Wkf{xi;6k), may be in- 
terpreted two ways. We may consider it to be the value of a single, complex 
density function that happens to representable as a weighted sum. Alterna- 
tively, we may consider it to represent a choice of one of the K components, 
with probability Wk, as the source for the datum, which is then drawn from 
the component density /(xj; Ok)- We can make the latter interpretation more 
explicit by introducing latent labels, {A,}, with Xi = k denoting assignment 
of datum i to component k. To keep track of the probabilities for particular 
assignments in equation (29), we use the labels to distinguish the values 
of Wk associated with the various datum factors, by replacing Wk in the ith 
factor with wx. (and summing over values of Aj rather than k) . Abbreviating 
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fixi] Ok) by fk,i-, we can rewrite the likelihood as follows: 



(30) C{{wk},{Ok}) = \{ 




where for the last hne we have collected factors of a particular weight by 
introducing a multiplicity function, mfc(A), counting the number of times 
component index k appears in the list of A'^ labels Aj. This dual representa- 
tion is well known in the literature on FMMs (see, e.g., [6]). 

Turning now to the cosmic ray likelihood function in equation (19), the 
product factor has the algebraic form of a FMM likelihood, with the source 
fluxes, Ffc, playing the role of mixing weights (but now no longer normal- 
ized). Equation (20) then follows by the same manipulations as shown above, 
with the addition of splitting the exponentiated sum in equation (19) into 

separate e"^'-'^'' factors, grouped with with their associated factors. 
The resulting likelihood function essentially corresponds to a Poissonized 
FMM. 

APPENDIX C: COMPUTATIONAL METHODS 

C.l. Algorithm for Markov chain Monte Carlo. To draw poste- 
rior samples, we perform Metropolis-within-Gibbs sampling on parameters 
/, Ft and A, using Gibbs sampling for Ft and A, and Metropolis sampling 
for /. The Gibbs sampling steps alternate between sampling from the full 
conditional distribution for Ft (i.e., the distribution given the data and all 
other parameters), and that for A. The full conditionals may be derived from 
(21) and the (Ft,/) priors. The conditional for the total flux is 



(31) Ft!/, a, L> ~ Gamma \ Nq + I 



1 



' 1 



+ (l-/)eo + /Efc>i"^fcefc 



where Gamma(a, s) denotes the gamma distribution with shape parameter 
a and scale parameter s. (Note that this distribution happens to be inde- 
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pendent of A.) The conditional for the cosmic ray labels is a multinomial 
distribution with probabilities 

(32) P{Xi\FTj,D)cx^xhx^, where hj = r}~^^'' V-Z^' ■ 

ex, [jWjej if J > 1. 

Finally, the conditional for the associated fraction is 



P{f\X,FT,D) a exp<^-FT 



(l-/)eo + /V 



k>l 



ekWk 



^33^ x(l — fy'^oW+b-l jNc-mo{X)+a-l 

Ft and A are sampled directly from the gamma and multinomial distribu- 
tions. / is sampled using a random walk Metropolis algorithm with Gaussian 
proposals centered around the current value of /. The variance of the Gaus- 
sian proposal density was tuned so that the acceptance rate was about 25%. 

C.2. Marginal Likelihood and Bayes Factor Computation. Fol- 
lowing Chib (1995), we can write the marginal likelihood for k as 

(34) ^^^^^^P{D\F^f*,X*)P{y\F^J*)P{F^)P{n 



P{F*J*,X*\D) 

where the double solidus indicates all the probabilities additionally condi- 
tion on K. Here F^,f*,X* are in principle arbitrary, but in practice should 
correspond to a point with high posterior density. All the terms in the nu- 
merator can be computed analytically, using the priors and the likelihood 
from equation (22). The denominator can be expressed as 

(35) P{F^,f*,X*\D) = P{f*\F},X*,D)P{F^\X*,D)P{X*\D) 

The first term on the right hand side is simply the full condition of /* 
evaluated at F^ and A*. Note that the normalizing constant can be computed 
using numerical integration. The remaining two terms need to be estimated 
using MCMC and can be done as follows: 

1 ^ 

(36) P(F^ |A*, K, D)^-Y^ P{F^\f'^^\X*,^, D), 

3=1 

1 ^ 

(37) P{X*\^,D) ^-J2P{^*\f^'\4'\^^D). 

9=1 
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Here G denotes the number of iterations. /^^^ and F!f denote the sample 
from the MCMC in iteration g. f'^^^ is a sample from a new MCMC run 
using the full conditionals given earlier with A fixed at A* in iteration g. 

For each k of interest, we first ran 5,000 iterations of Metropolis within 
Gibbs to obtain the high posterior density values, F^, /*, A*. Then, for each 
K, 3 additional chains of Gibbs sampling were run, each with 10,000 iter- 
ations. For subsequent calculations and plots, the chains were thinned, so 
that the lag-one autocorrelation is at most 0.15 for all parameters FT,f,X. 

We diagnosed convergence by visually inspecting the trace plots from dif- 
ferent chains, and by computing the Gelman- Rubin potential scale reduction 
statistic for chains of samples of continuous parameters such as / and Ft- 
For the discrete A parameters, we computed the fraction of time that each 
cosmic ray is assigned to each source during every 10 iterations. The diag- 
nostics were done on these fractions similarly to the continuous parameters. 

To validate our algorithms (including our convergence criteria) we devel- 
oped an enumerative algorithm that can directly calculate several poste- 
rior quantities for simplified models of small catalogs in the small-deflection 
regime, via a guided traversal of a tree of possible associations that has been 
thresholded to eliminate associations with negligible probability. We com- 
pared results from this deterministic algorithm with our MCMC results. We 
also used simulated datasets to verify that marginal distributions produce 
credible intervals of probability P that have prior-averaged coverage equal 
to P, a simplified version of the validation tests proposed by Cook, Gelman, 
and Rubin [14]. 

APPENDIX D: GEN A SINGLE-SOURCE MODEL 

Recent PAO results on the chemical composition of UHECRs, cited above, 
suggest that UHECRs may be predominantly heavy nuclei, which would suf- 
fer large magnetic defiections. Some investigators have suggested that UHE- 
CRs are all heavy nuclei from a single source — the nearest AGN, Cen A — 
with the apparent approximate isotropy of arrival directions a consequence 
of strong deflection [7, 17, 8]. 

As a simple test of this idea, we studied a model corresponding to our 
buckshot AGN association model, but with a single candidate source, Cen A, 
and with the association fraction / = 1. Figure 11 shows the Bayes factor 
comparing such a model to the isotropic background model, conditional on 
K, for various partitions of the data. For k = 0, the Cen A model makes the 
same predictions as the background model, and the Bayes factor is unity. As 
K increases, the Bayes factor never grows significantly larger than unity, so 
the Cen A model is not supported by any partition of the data, for any k. 
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Fig 11. Bayes factors comparing a model Cen A to be the source of all UHECRs agamst 
the null isotropic background model, conditional on k, shown as a function of k. Results 
are shown for various partions of the data (identified by line style, identified in the legend). 
Right panel expands the Iow-k region. 



For partitions of the data including the 42 Period 3 events, the Bayes factor 
opposes the Cen A model, strongly ruling out models with k>0.5, i.e., with 
deflection angular scales < 90°. 

The 90° scale is consistent with the expected angular scale for deflection 
by a regular magnetic field in equation (7) as long as Z is large (^15); but 
note that this equation is valid only in the small-deflection limit. The 90° 
scale is consistent with turbulent magnetic deflection only for the heaviest 
nuclei (e.g., for iron Z = 26) and for large magnetic field and length scales. 
For these astrophysical parameters, the Cen A model remains disfavored, 
although not overwhelmingly. We thus consider the PAO-10 data to argue 
against the hypothesis that all UHECRs originate from Cen A. 

APPENDIX E: COMPARISON WITH PRIOR BAYESIAN WORK 

Watson, Mortlock, and Jaffe [38] (WMJll) have independently devel- 
oped an approximate Bayesian approach for assessing evidence for associa- 
tion of the PAO-08 data (with data for the 27 UHECRs in periods 1 and 
2) with nearby AGN. They consider the best-fit CR directions as points 
drawn from a Poisson intensity on the celestial sphere constructed as a mix- 
ture of distributions located at AGN directions, with weights refiecting the 
AGN distances (corresponding to a standard-candle source intensity), and 
an isotropic background component. The source components have the form 
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of a "two-dimensional Gaussian on the sphere," with the probabihty density 
for a best-fit measured direction h arising from an AGN at vj given by 



with a dubbed the "smearing angle." This is a Fisher distribution (as in 
equation 9), with concentration parameter k = l/o"^. They focus on a model 
with o" = 3° (~ 0.052 rad) [sic], intended to reflect a combination of ~ 1° 
measurement uncertainties and few-degree magnetic deflections. They also 
provide a few results for a = 6° and 10° [sic].^'^ They calculate an approxi- 
mate likelihood function by finely pixelizing the celestial sphere, calculating 
the number of cosmic ray direction measurements expected in each pixel, 
and multiplying Poisson counting probabilities for the bins (with one count 
for bins containing a best-fit cosmic ray direction, and zero counts for the 
remaining bins). 

WMJll adopt a similar candidate host population as was used in the 
PAO-07 and PAO-08 analyses (the nearby AGN in the 12th VCV compi- 
lation; WMJll consider ~ 900 AGN within 100 Mpc) on the presumption 
that it is almost complete for the nearest AGN. They conclude that there 
is "strong evidence of a UHECR signal from the known VCV AGNs," such 
that at least some UHECRs come from AGN in the VCV catalog (or from 
sources within a few degrees of the AGN). For a 3° smearing angle, the 
marginal posterior density for the fraction associated with AGN (our / pa- 
rameter) has a mode of 0.15, and 68% highest density credible interval of 
[0.08,0.25]. For 6° and 10° smearing angles the estimated association frac- 
tion is larger, but the marginal distribution is also somewhat broader. They 
do not compare their association model to an alternative, and thus do not 
compute Bayes factors. 

Our analysis framework and our results differ in significant ways from 
those of WMJll (focusing on our results for the data from Periods 1 and 2). 
Methodologically, our approach is based on explicit modeling of associations 
(via marginal likelihood factors associating a particular cosmic ray with a 
particular AGN or the background), rather than considering best-fit cosmic 
ray directions to be samples from a point process. As noted in § 3.5, a factor 
in the likelihood function in our approach is analogous to that underlying 

^^We note that the parameter a is not an angle and does not have angular units 
(degrees or radians). We presume these dimensionally inconsistent equations imply solving 
the nonlinear equation relating a (or k) to the stated angular scale; e.g., for a 68.3% 
confidence or credible region with angular radius 6, in the small-angle limit a ~ 0.666 
(e.g., a ~ 0.035 for 3° uncertainties; see equation (11)). 



(38) 
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FMMs, but this mixture-like factor arises as a consequence of some of our 
modeling assumptions; it is not a starting point, and it does not hold for all 
astrophysically interesting association models that our framework accom- 
modates. For example, in § 6 we describe a more realistic family of magnetic 
deflection models ("radiant" models with exchangeable rather than IID de- 
flections for cosmic rays comprising a multiplet) whose likelihood function 
does not have the simple FMM form. In other common astronomical co- 
incidence assessment problems, such as establishing associations of sources 
detected in different electromagnetic wavebands, only singlet associations 
are meaningful; such models have no FMM representation but may be ac- 
commodated by our framework. Further discussion of this is in [28]. 

Another point of departure in methodology is that we distinguish mea- 
surement error from magnetic deflection. In the buckshot model adopted 
here, both the measurement error and deflection distributions are Fisher dis- 
tributions (note that the composition of these distributions is not a Fisher 
distribution). In radiant models, for example, the deflection distribution is 
more complicated. Explicit, separate treatment of these physically distinct 
effects enables handling heteroskedastic measurement errors for the UHECR 
directions (PAO-10 reports only a typical measurement error, but future 
catalogs will hopefully report the heteroskedastic uncertainties found in de- 
tailed air shower fits). Heteroskedastic uncertainties further thwart a simple 
FMM representation for the likelihood, again emphasizing the need for a 
framework built on explicit modeling of individual associations. 

Our approach also provides explicit estimates of probabilities for possible 
associations. In our MCMC algorithm, these can be found by calculating 
frequencies for different values of the A labels; we report such results in 
§ 4. WMJll report a weight for candidate associations, but note it is not a 
rigorous probability. 

WMJll analyze the data from Periods 1 and 2 jointly, without comment- 
ing on the possible effects of tuning on the implications of the Period 1 
data. 

Turning to astrophysical differences, we adopt the GIO [18] volume-complete 
catalog of nearby AGN as a candidate host catalog, rather than the 12th 
VCV catalog used by WMJll and others. Notably, 6 of the 17 AGN in the 
GIO catalog are not in the 12th VCV catalog (one of them appears in the 
more recent 13th VCV catalog). We chose the GIO catalog both for its com- 
pleteness, and because use of a small catalog was convenient for an initial 
study, as it enabled more extensive analyses of real and simulated data than 
would be possible with AGN from VCV. Figure 12 shows our marginal pos- 
terior density for / for k = 1000 (corresponding to a deflection scale ~ 2.7°) 
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Fig 12. Posterior distributions for f for model Mi, conditioned on ft = 1000, for various 
subsets of the PAO-10 data. 

for model Mi, for different subsets of the PAO-10 data. The mode based on 
the combined data from Periods 1 and 2 is at / ~ 0.1, about a two-thirds 
of the WMJll value of ~ 0.15, even though their catalog contains more 
than 50 times as many potential counterparts. Our common assumption of 
a standard candle intensity distribution is probably the main reason that 
the results are not more discrepant. In particular, although WMJll include 
AGN at distances to 100 Mpc in their catalog, the standard candle assump- 
tion forces the analysis to assign negligible detectable rates to all but the 
closest few AGN. The similarity of our estimates despite the disparity be- 
tween our AGN catalog sizes highlights how restrictive the standard candle 
assumption is. 

We explore a far greater range of magnetic deflection angular scales than 
did WMJll. Their discussion of deflection scales implicitly presumes UHE- 
CRs are light nuclei. Recent cosmic ray data and theoretical models moti- 
vate serious consideration of the possibility that many or most UHECRs are 
heavy nuclei, as noted above. 

Finally, we differ qualitatively in our conclusions about the strength of 
evidence for association of UHECRs with nearby AGN, particularly after ex- 
amining period-to-period differences, and considering the impact of period 3 
data (unavailable to WMJll). We quantify the support for association hy- 
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potheses by calculating Bayes factors explicitly comparing association and 
null models, both conditional on k (in Figure 4) and marginalized over a 
broad k range (in Table 2). 

WMJll do not calculate Bayes factors comparing their association and 
null models. Their claim of strong evidence for association appears to be 
based on the small marginal posterior density for values of the association 
fraction near zero. But this fails to distinguish parameter estimation from 
model assessment. As long as the likelihood for models with / = does not 
vanish, parameter estimation (with a continuous prior on /) simply does 
not address how the / = hypothesis compares to alternatives. To do this 
requires assigning a finite prior probability for / = 0, which we do here 
by considering it as a separate model and calculating Bayes factors. The 
Bayes factor comparing nested models depends on the size of the parameter 
space of the larger model in a way that accounts for "fine tuning" of the 
additional model parameters: the larger model will have parameter values 
producing better fits than the smaller model, but if the values of the addi- 
tional parameters are close enough to the default values corresponding to the 
smaller model, the marginal likelihood for the larger model will be smaller 
than that for the smaller model (the well-known "Ockham's razor" behav- 
ior of Bayes factors). Focusing on low-dimensional marginal distributions, 
such as the posterior density for /, can give an exaggerated impression of 
the strength of evidence for the larger model because it suppresses the large 
volume of parameter space associated with its additional parameters. Here, 
association models have not only the / parameter, but also many latent la- 
bel parameters (i.e., many association hypotheses that cannot be ruled out 
a priori). Calculation of the Bayes factor takes all of this into account. Using 
the Period 1 and Period 2 data available to WMJll does in fact produce 
large Bayes factors favoring association. But partitions of the data exclud- 
ing Period 1 data produce much smaller Bayes factors, even though the p{f) 
distributions found with these partitions assign very small density to / = 0. 
The Bayes factor calculations indicate that the complexity of association 
models may not be justified by existing data. 

Perhaps most importantly from an astrophysical perspective, we per- 
formed more extensive checking of our models, calling into question the 
adequacy of our shared isotropic background and standard candle assump- 
tions. We discuss this further in § 6. 

APPENDIX F: CHANGE-POINT ANALYSIS 

We can use the marginal likelihoods for the various periods to investigate 
whether parameter values for any of our models vary between time periods. 
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Fig 13. Bayes factors for change-point models. 



We compare models that allow the parameters to change between periods to 
models that keep the parameters the same for all periods. That is, we explore 
change-point models, with the change point locations at period boundaries. 
Specifically, we compute 3 quantities: 



^(l)(23) 



23 



B 



(2)(3) 



(39) 



where C 



(1)(2)(3) 



-^123 
-^23 



123 



is the marginal likelihood computed using the Chib estimate 
based on data from periods ii, . . . ,ig. For example, -B(i)(23) compares a model 
allowing parameters to differ between period 1 and the later periods, to a 
model with common parameter values across all periods. Note that -B(2){3) 
considers only untuned data. We compute -B(i)(23), i?(2)(3) and i?(i)(2){3) for 
both the 17 AGN and 2 AGN association models. Figure 13 shows these 
Bayes factors as functions of k. 

In the case of 17 AGN, for change-point models considering all of the data, 
the Bayes factors stay within [1/3,3] indicating no preference for one model 
over the other. The same is true for the 2 AGN model, except for k, > 300, 
where there is a modest preference for models with consistent parameter 
values across all periods. 

We also considered change point models that partition the data between 
(joint) association models and the null model, aiming to assess the possi- 
bility that the data are consistent with isotropy in some intervals, but with 
association (possibly spurious) in others. We again found Bayes factors to 
be equivocal. 
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